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Abstract 

The need for accurate material models to simulate the deformation, damage and failure of 
polymer matrix composites is becoming critical as these materials are gaining increased usage 
in the aerospace and automotive industries. While there are several composite material models 
currently available within LS-DYNA® , there are several features that have been identified that 
could improve the predictive capability of a composite model. To address these needs, a 
combined plasticity and damage model suitable for use with both solid and shell elements is 
being developed and is being implemented into LS-DYNA as MAT_213. A key feature of the 
improved material model is the use of tabulated stress-strain data in a variety of coordinate 
directions to fully define the stress-strain response of the material. To date, the model 
development efforts have focused on creating the plasticity portion of the model. The Tsai-Wu 
composite failure model has been generalized and extended to a strain-hardening based 
orthotropic material model with a non-associative flow rule. The coefficients of the yield 
function, and the stresses to be used in both the yield function and the flow rule, are computed 
based on the input stress-strain curves using the effective plastic strain as the tracking variable. 
The coefficients in the flow rule are computed based on the obtained stress-strain data. The 
developed material model is statable for implementation within LS-DYNA for use in analyzing 
the nonlinear response of polymer composites. 


Introduction 

As composite materials are gaining increasing use in aircraft components where impact 
resistance under high energy impact conditions is important (such as the turbine engine fan case), 
the need for accurate material models to simulate the deformation, damage and failure response 
of polymer matrix composites under impact conditions is becoming more critical. Within LS- 
DYNA 1 ' [1], several material models are available for application to the analysis of composites. 
For example, the Chang-Chang failure model [2] is utilized in MAT 22 and MAT 54. In these 
models, combinations of ratios of stresses to failure strengths are utilized to predict fiber or 
matrix based failure. The response is assumed be linear elastic, with limited capability to 
simulate the nonlinear shear response. In MAT 22 the failure is assumed to be brittle, while in 
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MAT 54 the composite elastic constants are selectively reduced based on the failure mode, and a 
gradual unloading is permitted until ultimate element failure is reached. In MAT58, a 
continuum damage model developed by Matzenmiller et al [3] is employed, where the initiation 
and accumulation of damage is assumed to be the primary driver of nonlinearity in the composite 
response. The failure stresses and strains of the material in each of the coordinate directions are 
specified by the user, and the material stress-strain curves are approximated based on this data. 
The original version of the material model assumes a material response that is independent of 
strain rate. However, in the enhanced model MAT158, a modified version of the model used in 
MAT 58 is employed which incorporates strain rate dependence into the material response [1]. 

A viscoelastic Prony series based on shear moduli is used to modify the computed stresses. The 
strain rate dependence in all of the coordinate directions is assumed to be identical. MAT161 is 
an additional continuum damage mechanics model available within LS-DYNA [4]. In this 
model, a stress to strength ratio approach in a variety of coordinate directions is used to specify 
the beginning of either fiber or matrix based failure. Appropriate failure mode based damage 
functions are then used to compute the reduction in elastic moduli in each of the coordinate 
directions. In MAT 219, based on the CODAM model developed at the University of British 
Columbia[5], a strain versus failure strain ratio based approach is used to predict the initiation of 
damage in a sublaminate of the composite. Separate damage accumulation and modulus 
reduction functions based on failure mode and coordinate direction are implemented in a 
continuum damage mechanics formulation. In MAT 22 1 , damage accumulation functions based 
on current strains, damage initiation strains and failure strains are used to reduce the elastic 
moduli of the composite in each of the coordinate directions in a damage mechanics approach 
[1]. In MAT 261, which is based on a damage and failure model developed by Pinho et al [6,7], 
separate models for fiber tension failure, fiber kinking failure, matrix tensile failure and matrix 
compressive failure are developed, where different functional forms are used for each of the 
failure criteria (for example, the matrix compression failure criterion is based on an extension of 
the Mohr-Coulomb failure criterion). The various failure criteria are combined together using 
fracture mechanics concepts to develop a constitutive model. In MAT 262, based on a model 
developed by Camonho, et al [8,9], an energy approach is utilized to generate damage functions 
in various coordinate directions within the context of a continuum damage mechanics 
formulation. 

While all of the material models discussed above have been utilized with some level of success 
in modeling the impact response of polymer composites, there are some areas where the 
predictive capability can be improved. In general, the existing models often require significant a 
priori knowledge of the damage and failure such that their use as predictive tools can be limited. 
While these models generally assume that the nonlinear response of the composite is due to 
either deformation mechanisms (plasticity) or damage mechanisms, an improved model should 
have the capability to simulate the actual material behavior in which the material nonlinearity is 
due to a combination of both deformation and damage mechanisms. The input to the current 
material models generally consists of point-wise properties that lead to curve fit approximations 
to the material stress-strain curves. This type of approach leads to either models with only a few 
parameters, which provide a crude approximation at best to the actual stress-strain curve, or to 
models with many parameters which require a large number of complex tests to characterize. An 
improved approach would be to use tabulated data from a well-defined set of experiments to 
accurately define the complete stress-strain response of the material. Furthermore, many of the 
existing models are only suitable for use with two-dimensional shell elements, which cannot 
capture the through- thickness response which may be significant in impact applications. Ideally, 
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a fully three-dimensional formulation suitable for use with solid elements would be available, 
along with the shell element formulation. In addition, incorporating the effects of strain rate is 
only possible in a few of the existing models, and in those models, the effects of strain rate on 
the material response are assumed to be the same in all of the material coordinate directions. 

To begin to address these needs, a multi-institution consortium has been formed to develop and 
implement a new composite material model within LS-DYNA, which will be implemented as 
MAT213. Currently, the primary focus of the effort has been on the development and testing of 
the orthotropic macroscopic three-dimensional plasticity model. As will be discussed in more 
detail in the following sections of this paper, the commonly used Tsai-Wu composite failure 
criterion [ 1 0] has been generalized and extended to a strain-hardening model with a non- 
associative flow rule. A similar approach has been used by Sun and Chen [1 1], in which a 
general quadratic plastic potential function was generalized into a plasticity model suitable for 
use with composites. The coefficients of the yield function for the new composite model are 
determined based on tabulated stress-strain curves in the various normal and shear directions, 
along with selected off-axis curves. The non-associative flow rule is used to compute the 
components of the plastic strain along with the effective plastic strain. The stresses that are used 
in the flow rule are computed based on the input tabulated stress-strain curves. Systematic 
procedures have been developed to compute the various coefficients in the yield function and the 
flow rule based on the tabulated input data. An important point to note is that the developed 
material law is not limited to the analysis of unidirectional plies alone. Alternatively, the 
material model is meant to be a fully generalized model suitable for use with any composite 
architecture (unidirectional, laminated or textile). In the following sections of this paper, a 
detailed derivation of the yield function and the flow law are presented. The procedures to be 
used to characterize the material constants in the yield function and the flow law are also 
discussed in detail. Details of the numerical implementation of the material model within LS- 
DYNA and descriptions of verification and validation studies that were conducted using the 
developed material model are available in a companion paper [12]. 

Material Law Derivation 


A general quadratic three-dimensional ortho tropic yield function based on the Tsai-Wu failure 
model is specified as follows, where 1, 2, and 3 refer to the principal material directions. 
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In the yield function, cry represents the stresses and F l} and F k are coefficients that vary based on 
the current values of the yield stresses in the various coordinate directions. By allowing the 
coefficients to vary, the yield surface evolution and hardening in each of the material directions 
can be precisely defined. The values of the normal and shear coefficients can be determined by 
simplifying the yield function for the case of unidirectional tensile and compressive loading in 
each of the coordinate directions along with shear tests in each of the shear directions, with 
results as shown below. 
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In the above equation, the stresses are the current value of the yield stresses in the normal and 
shear directions (detennined using procedures to be discussed below), where the superscript T 
indicated the tensile yield stress, and the superscript C indicates the absolute value of the 
compressive yield stress. To determine the values of the off-axis coefficients (which are required 
to capture the stress interaction effects), the results from 45° off-axis tests in the various 
coordinate directions can be used. For example, the stresses in the local material axis system 
resulting from a uniaxial tensile test of a [45°] composite (for a unidirectional composite), or any 
uniaxial tensile test conducted at 45° in the 1-2 plane from the longitudinal (1) material axis for a 
multi-ply laminated or textile composite, can be computed to be the following by the stress 
transformation equations [10]. 

ffii = 0.5crj, 45 

(T 2 2 = 0-5cr y45 (3) 

(J l2 = 0.5cr r45 


In the above equation, oy 45 is the structural level current uniaxial yield stress resulting from a 
tensile test on a [45°] ply, and the remaining stresses are the material axis system based stresses. 
By substituting this equation into Equation (1) and solving for the off-axis constant Fn, the 
following expression can be obtained 
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where all of the terms are as defined earlier. By using similar procedures involving 45° rotations 
in the 2-3 plane and the 1-3 plane, similar expressions for the constants F 23 andfffi can be 
determined. 
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A non-associative flow rule is used to compute the evolution of the components of plastic strain. 
The plastic potential for the flow rule is shown below 


h — Ill'll + H 22 CT 21 + ^ 33^33 + 2.H l2 CT\ j<T 22 + 2// 23 C7 22 C7 3 3 + 2H 3l O’ 33 (J n +H 44 CT l2 + // ? 5 <T 2 3 +H 66 ct 31 (5) 


where oy are the current values of the stresses and //,, are independent coefficients, which are 
assumed to remain constant. Procedures for determining the values of the coefficients will be 
discussed later in this paper. The plastic potential function in Equation (5) is used as follows in a 
flow law, where the usual normality hypothesis from classical plasticity [13] is assumed to apply 
and the variable X is a scalar plastic multiplier. 
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In the above equation, s p are the components of plastic strain rate. For the shear components, 

tensorial, not engineering, strains are used. By examining ratios of the components of plastic 
strain in response to a given uniaxial stress component, the normal and off-axis coefficients H l} in 
the plastic potential function can be related to “plastic Poisson’s ratios”, Vy P , as follows, which as 
will be shown later in this paper, is useful in developing procedures to characterize the values of 
the coefficients. 
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In the above equation the expressions in the first column are only valid for the case where the 
only non-zero stress is in the 1 direction, the expressions in the second column are only valid for 
the case where the only non-zero stress is in the 2 direction, and the expressions in the third 
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column are only valid for the case where the only non-zero stress in in the 3 direction. An 
important point to note is that while the yield function can accommodate differences between the 
tensile and compressive responses of the composite, the flow law, due to its lack of linear terms, 
cannot make such a distinction. However, introduction of a linear term would make the plastic 
Poisson’s ratio, and thus the coefficients in the flow law, highly dependent on the stress which 
would most likely lead to erratic behavior of the model. Another consideration is that the 
coefficients of the flow law could be assumed to vary based on the current stress and strain state, 
but that would require developing evolution conditions on the flow law, which would be very 
difficult to detennine based on experimental data. 

Given the flow law, the principal of the equivalence of plastic work [13] can be used to 
detennine expressions for the effective stress and effective plastic strain. By taking the vector 
product of the stress tensor and the plastic strain rate tensor, one ends up with the plastic 
potential function multiplied by the plastic multiplier A . Since by the principal of the 
equivalence of plastic work the vector product of the stress tensor and the plastic strain rate 
tensor must equal the product of the effective stress and the effective plastic strain rate, one can 
conclude that the plastic potential function h can be defined as the effective stress and the plastic 
multiplier X can be defined as the effective plastic strain rate. This concept is expressed 
mathematically in the equation shown below 


where a e is the effective stress and s p e is the effective plastic strain rate. 

To compute the current value of the yield stresses needed for the yield function, the common 
practice in plasticity constitutive equations is to use analytical functions to define the evolution 
of the stresses as a function of the components of plastic strain (or the effective plastic strain). 
Alternatively, in the developed model tabulated stress-strain curves are used to track the yield 
stress evolution. The user is required to input twelve stress versus plastic strain curves. 
Specifically, the required curves include uniaxial tension curves in each of the normal directions 
(1,2,3), uniaxial compression curves in each of the nonnal directions (1,2,3), shear stress-strain 
curves in each of the shear directions (1-2, 2-3 and 3-1), and 45 degree off-axis tension curves in 
each of the 1-2, 2-3 and 3-1 planes. The 45 degree curves are required in order to properly 
capture the stress interaction effects. By utilizing tabulated stress-strain curves to track the 
evolution of the deformation response, the experimental stress-strain response of the material can 
be captured exactly without any curve fit approximations. 

The required stress-strain data can be obtained either from actual experimental test results or by 
appropriate numerical experiments utilizing stand-alone codes. Currently, only static test data is 
considered. Future efforts will involve adding strain rate and temperature dependent effects to 
the computations. For certain composite architectures, the number of stress-strain curves that 
need to be input can be simplified. For example, for unidirectional composites, the normal 
response in the 3 direction can be assumed to be identical to the normal response in the 2 
direction due to the transverse isotropy of the material. Likewise, the shear response in the 1-3 
direction can be assumed to be identical to the shear response in the 1-2 direction. Similarly, the 
45 degree off-axis curve in the 1-3 plane can be assumed to be identical to the 45 degree off-axis 
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curve in the 1-2 plane. Due to the transverse isotropy of a unidirectional composite, one can also 
demonstrate that for the case of a unidirectional composite the 45 degree off-axis curve in the 2-3 
plane is approximately equal to the normal stress-strain curve in the 2 (or 3 direction). For other 
composite architectures with significant symmetry (such as two-dimensional plain or satin weave 
composites), a certain degree of simplification in the stress-strain data that is required to be input 
is also possible. 

To track the evolution of the deformation response along each of the stress-strain curves, the 
effective plastic strain is chosen to be the tracking parameter. Using the numerical procedure 
described briefly later in this paper and in more detail in a companion paper [12], the effective 
plastic strain is computed for each time/load step. The stresses for each of the tabulated input 
curves corresponding to the current value of the effective plastic strain are then used to compute 
the yield function coefficients. To enable this tracking, each of the input stress-strain curves 
needs to be converted into a plot of stress versus effective plastic strain. To enable this 
conversion, the principal of the equivalence of plastic work is used again to relate the 
unidirectional plastic strain increment to the effective plastic strain increment. For example, for 
the case of unidirectional loading in the 1 direction the effective plastic strain increment (and 
integrated current total plastic strain) can be computed as follows 
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where cq i is the unidirectional stress in the loading direction, ds[\ is the plastic strain increment 
in the loading direction and ds % is the increment in effective plastic strain. Note that in a 
general case the increment in effective plastic strain can be computed as the incremental area 
under the stress-plastic strain curve divided by the current value of the effective stress. A 
graphical representation of process is shown below in Figure 1 . 
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User provided load curves are Internally stored is true stress 

true stress versus true plastic strain versus effective plastic strain 

Figure 1 : Conversion of stress versus plastic strain curves to stress versus effective plastic strain 

curves. 


To compute the evolution of the effective plastic strain, a numerical algorithm based on the 
radial return method is employed. Given the revised value of the effective plastic strain, the 
yield stress values and the overall stress state of the material can be updated. Details of the 
numerical implementation are provided in a companion paper [12], but some of the key 
theoretical fundamentals of the algorithm are given below. The standard elastic constitutive 
equation is used to compute the revised stresses, where the flow law is applied for the 
computation of the plastic strains. 
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In the above equation, C is the standard elastic stiffness matrix, e is the total strain rate tensor, 
and the remaining terms are as defined previously. To compute the effective plastic strain rate, 
the consistency condition is applied in combination with the elastic constitutive equation 
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where /is the yield function defined in Equation (1) and q is the vector of yield stresses based on 
the input stress-strain curves. 
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where <T r45 „ 23 is the yield stress resulting from a uniaxial tensile test conducted at 45° in the 2-3 
plane and ov 45 _ 31 is the yield stress resulting from a uniaxial tensile test conducted at 45° in the 

3-1 plane. To compute the derivative of the components of the q vector with respect to the 
effective plastic strain X, a chain rule type of computation is carried out where the derivatives of 
the individual components of the vector (the respective yield stresses) with respect to the 
effective plastic strain are determined as follows. 


da,, da., ds f 
dX ds.j dX 


(13) 


The derivative of the components of the yield stress vector with respect to plastic strain is the 
instantaneous slope of the stress versus plastic strain curve from each of the 12 input stress-strain 
curves. The derivative of the plastic strain with respect to the plastic multiplier X can be 
computed from the flow law for each of the specialized unidirectional (or off-axis) cases 
represented by the input stress-strain curves. Note that for the off-axis cases the flow law must 
first be written in the material axis system and then the plastic strains must be transformed into 
the structural axis system in order to compute the proper derivative of the plastic strain with 
respect to the plastic multiplier. A secant iteration approach is used in combination with the 
radial return method to compute the specific value of the effective plastic strain rate [12]. 


Characterization of Flow Law Coefficients 

The coefficients in the flow rule need to be characterized based on data obtained from 
experimental stress-strain curves for the composite. Specifically, data from the 12 experiments 
leading to the input stress-strain curves must be utilized. The data can be obtained either from 
actual experiments or numerical experiments. For example, if one knows the mechanical 
properties and stress-strain curves for the fiber and matrix constituents, numerical analyses can 
be conducted using either high fidelity finite element analysis (where the fiber and matrix are 
modeled as distinct constituents) or analytical tools such as the NASA Glenn developed 
micromechanics code MAC/GMC [14] to obtain numerical approximations to the actual 
composite experimental data. 

First, the procedures used to obtain the required coefficients for the simplified case of a 
unidirectional carbon fiber based polymer matrix composite will be discussed. Afterward, a 
generalized procedure to obtain the required constants for a general composite will be presented. 
For a unidirectional carbon fiber composite, as discussed by Sun and Chen [1 1] a reasonable 
approximation to make is that due to the linear behavior of the carbon fiber the plastic strain in 
the fiber direction (the 1 direction) is equal to zero for all values of stress. By examining the 
second expression in Equation (6), one can conclude that for a general stress state the only way 
the plastic strain in the 1 direction can be zero is for the coefficients H n , Hn and //] 3 to all be 
identically zero. Next, since the transverse tension response of the composite can display a 
certain degree of nonlinearity, a reasonable approximation to make is that for the case of a 
unidirectional load in the 2 direction the effective stress is equal to the applied transverse tension 
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stress 022- By simplifying the plastic potential function g (shown in Equation (5)) for the case of 
a unidirectional transverse tension stress, and by setting this term (along with the transverse 
tensile stress) equal to the effective stress, one can show that the value of the coefficient H 22 
must equal one. Since a unidirectional composite is transversely isotropic, one can set the value 
of H 33 to also equal one. Given these values, using the definition of the plastic Poisson’s ratio 
v 23 shown in Equation (7), one can find that the coefficient H23 equals the negative of the plastic 
Poisson’s ratio v 23 ■ Since the flow surface coefficients are assumed to be constant, an average 
value of the plastic Poisson’s ratio must be determined from the experimental and/or numerical 
test data from a unidirectional transverse tension test. To determine the in-plane shear 
coefficient H44, utilizing a procedure similar to that developed by Sun and Chen [11] for their 
model, the plastic potential function g given in Equation (5) and the plastic strain definition 
given in Equation (6) can be simplified for the case of in-plane shear loading in the 1-2 direction, 
yielding the following expressions. 
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Starting from the experimental stress versus plastic strain plot for the case of in-plane shear 
loading, using the expressions in Equation (14) the value of H44 can be optimized to provide the 
closest possible match with the overall effective stress versus effective strain curve, selected for 
the case of a unidirectional carbon fiber composite to be the transverse tension stress versus 
plastic strain curve. Due to the transverse isotropy of the unidirectional composite, the value of 
the coefficient H ^ can be set equal to the value of the coefficient H44. To determine the value of 
the If 55 coefficient, a similar procedure can be carried out for the case of shear loading in the 2-3 
direction. Alternatively, due to the transverse isotropy of the composite, by computing the 
effective stress for the case of off-axis loading in the 2-3 plane, an isotropic type of relation can 
be determined for the coefficient H 55 

H, s = 2(l + v') (15) 

where u 2 3 P is the plastic Poisson’s ratio in the 2-3 direction, obtained from the material data. 

Similar procedures are used for the case of a general composite with no identified symmetries. 
An example of a general composite would be a triaxially braided composite which is modeled as 
a smeared, homogenous material. Another example would be a generalized laminated composite 
modeled as a smeared, homogeneous material. For the case of a general composite, a 
unidirectional tension test in the longitudinal (1) direction is chosen as the baseline case. By 
utilizing Equation (5) and setting the plastic potential function and the unidirectional longitudinal 
stress equal to the effective stress, the value of the H n coefficient can be found to equal one. 
Based on the relations given in Equation (7), the value of the coefficient Hn can be found to 
equal the negative of the plastic Poisson’s ratio v p 2 , and the value of the coefficient /ffi can be 
found to be equal to the negative of the plastic Poisson’s ratio vf 3 . As before, the value of the 
plastic Poisson’s ratio used in the computations is the average value of the parameter determined 
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from the experimental and/or numerically obtained data. In a similar fashion, the values of the 
coefficients H 2 2 , H 23 and H 33 can be computed. 


To compute the values of the coefficients H44, // 55 , and H 6(l , the same procedure that was utilized 
to compute H44 for the case of a unidirectional carbon fiber composite is used, only all three 
shear loading cases ( 012 , o 33 , and a 3 i) must be used. 


A generalized plasticity model suitable for use in polymer composite impact simulations has 
been developed. The plasticity model will be part of a generalized plasticity and damage model 
which will be implemented into LS-DYNA as MAT 2 13. The yield function is based on the 
Tsai-Wu composite failure model, and a suitable nonassociative flow rule was defined. Methods 
of utilizing tabulated stress-strain data to track the evolution of the yield stresses as a function of 
the effective plastic strain have been developed. The elastic constitutive equation can be used in 
combination with the consistency condition to compute the effective plastic strain. In a 
companion paper [ 12 ], further details of the numerical algorithm are provided along with 
examples of verification and validation studies which have been performed to examine the 
accuracy and efficiency of the analytical model. 

Future efforts will involve the development of an equivalent damage model, to be employed in 
combination with the plasticity model, which will also be based on utilizing tabulated stress- 
strain data to track the evolution of damage in the composite. Element failure will also be 
incorporated within the model. In addition, strain rate and temperature effects will be added to 
the analysis capability. Overall, when completed the composite model as implemented in 
MAT 213 will provide significant improvements to the state of the art in the modeling of the 
impact response of polymer matrix composites. 
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